home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / math / alged34.zip / ALGEDSRC.ZIP / ALGEXPR.C < prev    next >
C/C++ Source or Header  |  1996-06-06  |  20KB  |  787 lines
  1. /*--------------------------------------------------------------------
  2.    Alged:  Algebra Editor    henckel@vnet.ibm.com
  3.  
  4.    Copyright (c) 1994 John Henckel
  5.    Permission to use, copy, modify, distribute and sell this software
  6.    and its documentation for any purpose is hereby granted without fee,
  7.    provided that the above copyright notice appear in all copies.
  8. */
  9. #include "alged.h"
  10. /*--------------------------------------------------------------------
  11.    associate
  12.    rotate the association on add,sub,mul,div
  13.    a+b-c => a-c+b => b+a-c => b-c+a
  14.    note, we make use of the fact that ADD+1 = SUB and ADD is even.
  15. */
  16. void associate(node *p) {
  17.   node *a,*b;
  18.   int opr;
  19.  
  20.   opr = p->kind;
  21.   b = p->rt;
  22.  
  23.   if (opr==DIV && b->kind==MUL) {     /* special handle for bisected */
  24.     p->rt = b->rt;
  25.     b->rt = b->lf;
  26.     b->lf = p->lf;
  27.     p->lf = b;
  28.     b->kind = DIV;
  29.     b = p->rt;
  30.   }
  31.   if (opr==ADD || opr==MUL || opr==SUB || opr==DIV) {
  32.     a = p;
  33.     while ((a->lf->kind|1) == (opr|1)) {
  34.       a->kind = a->lf->kind;
  35.       a->rt = a->lf->rt;
  36.       a = a->lf;
  37.     }
  38.     a->kind = opr;
  39.     if (opr&1) {           /* subtr or divide */
  40.       a->rt = b;
  41.     }
  42.     else {
  43.       a->rt = a->lf;
  44.       a->lf = b;
  45.     }
  46.   }
  47.   else if (opr==EQU) {           /* commute equality */
  48.     p->rt = p->lf;
  49.     p->lf = b;
  50.   }
  51. }
  52.  
  53. /*--------------------------------------------------------------------
  54.    commute
  55.    commute on add,sub,mul,div
  56.    e.g.   a/b ==>  b^(-1)/a^(-1)
  57. */
  58. void commuteNOTUSED(node *p) {
  59.   node *a,*b;
  60.  
  61.   a = p->lf;
  62.   b = p->rt;
  63.   if (p->kind==MUL || p->kind==ADD || p->kind==EQU) {
  64.     p->lf=b; p->rt=a;
  65.   }
  66.   else if (p->kind==DIV || p->kind==SUB) {     /* non-commutative */
  67.     if (a->kind==p->kind+1 &&
  68.         b->kind==p->kind+1 &&
  69.         a->rt->kind==NUM &&
  70.         b->rt->kind==NUM) {              /* has coefficients */
  71.       p->lf=b; p->rt=a;
  72.       a->rt->value = -a->rt->value;
  73.       b->rt->value = -b->rt->value;
  74.     }
  75.     else {                          /* need to make coefficients */
  76.       p->lf = newoper(p->kind+1);
  77.       p->rt = newoper(p->kind+1);
  78.       p->lf->lf = b;
  79.       p->lf->rt = newnum(-1);
  80.       p->rt->lf = a;
  81.       p->rt->rt = newnum(-1);
  82.     }
  83.   }
  84. }
  85.  
  86. /*--------------------------------------------------------------------
  87.       a^x * a^y ==> a^(x+y)
  88.       x^a * y^a ==> (xy)^a
  89. */
  90. int expjoin(node *p) {
  91.   int i,r=0;
  92.   node *a,*b;
  93.  
  94.   for (i=0; i<p->nump; ++i)
  95.     r+=expjoin(p->parm[i]);
  96.  
  97.   a = p->lf;
  98.   b = p->rt;
  99.   if (p->kind==MUL || p->kind==DIV) {
  100.     if (a->kind==EXP && b->kind==EXP &&
  101.         equal(a->lf,b->lf)) {      /*  a^x * a^y = a^(x+y) */
  102.       freetree(a->lf);
  103.       a = a->rt;
  104.       freenode(p->lf);
  105.       b->kind = p->kind-2;
  106.       p->kind = EXP;
  107.       p->lf = b->lf;
  108.       b->lf = a;
  109.     }
  110.     else if (b->kind==EXP &&
  111.         equal(a,b->lf)) {               /*  a * a^y = a^(1+y) */
  112.       freetree(a);
  113.       a = newnum(1);
  114.       b->kind = p->kind-2;
  115.       p->kind = EXP;
  116.       p->lf = b->lf;
  117.       b->lf = a;
  118.     }
  119.     else if (a->kind==EXP &&
  120.              equal(a->lf,b)) {          /*  a^x * a = a^(x+1) */
  121.       freetree(b);
  122.       b = newnum(1);
  123.       a->kind = p->kind-2;
  124.       p->kind = EXP;
  125.       p->lf = a->lf;
  126.       a->lf = a->rt;
  127.       a->rt = b;
  128.       p->rt = a;
  129.     }
  130.     else if (a->kind==EXP && b->kind==EXP &&
  131.         equal(a->rt,b->rt)) {      /*  x^a * y^a = (xy)^a */
  132.       freetree(b->rt);
  133.       p->rt = a->rt;
  134.       a->rt = b->lf;
  135.       freenode(b);
  136.       a->kind = p->kind;
  137.       p->kind = EXP;
  138.     }
  139.     else return r;
  140.     ++r;
  141.   }
  142.   return r;
  143. }
  144.  
  145. /*--------------------------------------------------------------------
  146.    remove sub and div
  147.  
  148.    x/y  =  x*y^-1
  149.    x-y  =  x+y*-1
  150. */
  151. int nosubdiv(node *p) {
  152.   int i,r=0;
  153.   node *a,*b;
  154.  
  155.   for (i=0; i<p->nump; ++i)
  156.     r+=nosubdiv(p->parm[i]);
  157.  
  158.   a = p->lf;
  159.   b = p->rt;
  160.   if (p->kind==SUB || p->kind==DIV) {
  161.     ++r;
  162.     --p->kind;
  163.     if (b->kind==NUM && p->kind==ADD)           /* a-2 = a+(-2) */
  164.        b->value = -b->value;
  165.     else if (b->kind==NUM && !whole(b->value))
  166.        b->value = 1.0/b->value;
  167.     else {
  168.       a = newoper(p->kind+2);    /* MUL or EXP */
  169.       a->lf = b;
  170.       a->rt = newnum(-1);
  171.       p->rt = a;
  172.     }
  173.   }
  174.   return r;
  175. }
  176.  
  177. /*--------------------------------------------------------------------
  178.    bisect node
  179.  
  180.    This converts expressions to canonical form in which
  181.    1. + * are left assoc, ^ is right assoc.
  182.    2. scope of / is increased, e.g. (a/b)*c ==> ac/b,
  183.    3. x+y*-2 is changed to x-y*2
  184.    3' x+(-2) is changed to x-2
  185.    4. x*y^-2 is changed to x/y^2
  186.    4' x*(0.5) is changed to x/2
  187.    5. scope of ^ is reduced
  188.    The name "bisect" means that each MUL clag is broken into two pieces the
  189.    numerator elements and the denominator ones.
  190. */
  191. int bisect(node *p) {
  192.   int i,r=0;
  193.   node *a,*b;
  194.  
  195.   for (i=0; i<p->nump; ++i)
  196.     r+=bisect(p->parm[i]);
  197.  
  198.   a = p->lf;
  199.   b = p->rt;
  200.   switch (p->kind) {
  201.   case ADD:
  202.   case MUL:
  203.     /*--------------------------------------------------------------------
  204.        Add or Multiply
  205.     */
  206.     if (b->kind==p->kind ||              /* a+(b+c) = a+b+c */
  207.         b->kind==p->kind+1) {            /* a+(b-c) = a+b-c */
  208.       swingb;
  209.       i = b->kind;
  210.       b->kind = p->kind;
  211.       p->kind = i;
  212.     }
  213.     else if (p->kind==ADD &&             /* a+(-2) = a-2 */
  214.              b->kind==NUM &&
  215.              b->value < 0) {
  216.       p->kind = SUB;
  217.       b->value = -b->value;
  218.     }
  219.     else if (b->kind==p->kind+2 &&       /* a+b*-2 = a-b*2 */
  220.         b->rt->kind==NUM &&
  221.         b->rt->value < 0) {
  222.       ++p->kind;
  223.       b->rt->value = -b->rt->value;
  224.     }
  225.     /*--------------------------------------------------------------------
  226.        Add only
  227.     */
  228.     else if (p->kind==ADD &&
  229.         aop(b->kind)) {                   /* a+(b-c) = (a+b)-c */
  230.       swingb;
  231.       p->kind = b->kind;
  232.       b->kind = ADD;
  233.     }
  234.     /*--------------------------------------------------------------------
  235.        Multiply only
  236.     */
  237.     else if (p->kind==MUL &&
  238.         a->kind==p->kind+1) {       /* a-b+c = a+c-b */
  239.       p->rt = a->rt;
  240.       a->rt = b;
  241.       ++p->kind;
  242.       --a->kind;
  243.     }
  244.     else if (p->kind==MUL &&
  245.         a->kind==p->kind+2 &&       /* b*-2+a = a-b*2 */
  246.         a->rt->kind==NUM &&
  247.         a->rt->value < 0) {
  248.       ++p->kind;
  249.       p->lf = b;
  250.       p->rt = a;
  251.       a->rt->value = -a->rt->value;
  252.     }
  253.     else break;
  254.     ++r; break;
  255.   case SUB:
  256.   case DIV:
  257.     /*--------------------------------------------------------------------
  258.        Subtract or divide
  259.     */
  260.     if (b->kind==p->kind+1 &&       /* a-b*-2 = a+b*2   NOT NECES. */
  261.         b->rt->kind==NUM &&
  262.         b->rt->value < 0) {
  263.       --p->kind;
  264.       b->rt->value = -b->rt->value;
  265.     }
  266.     /*--------------------------------------------------------------------
  267.        Subtract only
  268.     */
  269.     else if (p->kind==SUB &&             /* a-(-2) = a+2 */
  270.         b->kind==NUM &&
  271.         b->value < 0) {
  272.       p->kind = ADD;
  273.       b->value = -b->value;
  274.     }
  275.     else if (p->kind==SUB &&
  276.          aop(b->kind)) {                   /* a-(b+c) = (a-b)-c */
  277.       swingb;
  278.       p->kind = (ADD+SUB) - b->kind;
  279.       b->kind = SUB;
  280.     }
  281.     /*--------------------------------------------------------------------
  282.        Divide only
  283.     */
  284.     else if (p->kind==DIV &&
  285.         b->kind==p->kind) {             /* a-(b-c) = a+c-b */
  286.       p->lf = b;
  287.       p->rt = b->lf;
  288.       b->lf = a;
  289.       --b->kind;
  290.     }
  291.     else if (p->kind==DIV &&
  292.         a->kind==p->kind) {        /* a-b-c = a-(b+c) */
  293.       swinga;
  294.       --a->kind;
  295.     }
  296.     else break;
  297.     ++r; break;
  298.   case EXP:
  299.     /*--------------------------------------------------------------------
  300.        exponent
  301.     */
  302.     if (a->kind==EXP) {                  /* (x^y)^z = x^(y*z) */
  303.       swinga;
  304.       a->kind = MUL;
  305.     }
  306.     else if (b->kind==NUM &&             /* a^(-2) = 1/a^2 */
  307.              b->value < 0) {
  308.       p->kind = DIV;
  309.       p->lf = newnum(1);
  310.       p->rt = newoper(EXP);
  311.       p->rt->lf = a;
  312.       p->rt->rt = b;
  313.       b->value = -b->value;
  314.     }
  315.     else break;
  316.     ++r; break;
  317.   }
  318.   return r;
  319. }
  320.  
  321. /*--------------------------------------------------------------------
  322.    exponent expand - expand any integer exponents less than 100.
  323. */
  324. int exexpand(node *p) {
  325.   int i,r=0;
  326.   node *a,*b;
  327.  
  328.   for (i=0; i<p->nump; ++i)
  329.     r+=exexpand(p->parm[i]);
  330.  
  331.   a = p->lf;
  332.   b = p->rt;
  333.   if (p->kind==EXP &&
  334.       b->kind==NUM &&           /* a^n = a^(n-1)*a */
  335.       b->value > 1 &&
  336.       b->value <= maxpow &&
  337.       whole(b->value)) {
  338.     if (--b->value == 1) {
  339.       p->lf = deepcopy(a);
  340.       freenode(b);
  341.     }
  342.     else {
  343.       b = newoper(EXP);
  344.       b->lf = deepcopy(a);
  345.       b->rt = p->rt;
  346.       p->lf = b;
  347.     }
  348.     p->rt = a;
  349.     p->kind = MUL;
  350.     ++r;
  351.   }
  352.   return r;
  353. }
  354.  
  355. /*-----------------------------------------------------------------
  356.    within - if p is a factor within q, then return the rest of q.
  357.    if the result is not null, then it is a new copy.
  358. */
  359. node *within(node *p, node *q) {
  360.   node *a=q,*b=NULL;
  361.  
  362.   while (a->kind==MUL) {
  363.     if (equal(p,a->rt)) {
  364.       if (b) {
  365.         b->lf = a->lf;             // skip over match
  366.         q = deepcopy(q);           // make copy
  367.         b->lf = a;                 // repair
  368.       }
  369.       else q = deepcopy(a->lf);
  370.       return q;
  371.     }
  372.     else if (equal(p,a->lf)) {
  373.       if (b) {
  374.         b->lf = a->rt;             // skip over match
  375.         q = deepcopy(q);           // make copy
  376.         b->lf = a;                 // repair
  377.       }
  378.       else q = deepcopy(a->rt);
  379.       return q;
  380.     }
  381.     b = a;
  382.     a = b->lf;
  383.   }
  384.   return NULL;
  385. }
  386.  
  387. /*--------------------------------------------------------------------
  388.    ComDeno - find common denominators.
  389. */
  390. int comdeno(node *p) {
  391.   int i,r=0;
  392.   double x;
  393.   node *a,*b,*nu,*de;
  394.  
  395.   for (i=0; i<p->nump; ++i)
  396.     r+=comdeno(p->parm[i]);
  397.  
  398.   a = p->lf;
  399.   b = p->rt;
  400.   if (aop(p->kind) && (a->kind==DIV || b->kind==DIV)) {
  401.     if (a->kind==DIV && b->kind==DIV) {
  402.       if (equal(a->rt,b->rt)) {     /* a/x + b/x = (a+b)/x */
  403.         freetree(b->rt);
  404.         de = a->rt;
  405.         a = a->lf;
  406.         b = b->lf;
  407.         freenode(p->lf);
  408.         freenode(p->rt);
  409.       }
  410.       else if (!!(nu=within(a->rt,b->rt))) {   /* a/x + b/xy = (ay+b)/xy */
  411.         freetree(a->rt);
  412.         a = cons(nu,MUL,a->lf);
  413.         de = b->rt;
  414.         b = b->lf;
  415.         freenode(p->lf);
  416.         freenode(p->rt);
  417.       }
  418.       else if (!!(nu=within(b->rt,a->rt))) {   /* a/xy + b/x = (a+by)/xy */
  419.         freetree(b->rt);
  420.         b = cons(nu,MUL,b->lf);
  421.         de = a->rt;
  422.         a = a->lf;
  423.         freenode(p->lf);
  424.         freenode(p->rt);
  425.       }
  426.       else {                        /* a/c + b/d = (ad+bc)/cd */
  427.         a->kind = b->kind = MUL;
  428.         de = b->rt;
  429.         b->rt = a->rt;
  430.         a->rt = de;
  431.         de = newoper(MUL);
  432.         de->lf = deepcopy(b->rt);
  433.         de->rt = deepcopy(a->rt);
  434.       }
  435.     }
  436.     else if (a->kind==DIV) {        /* a/b + c = (a + bc)/b */
  437.       a->kind = MUL;
  438.       de = a->lf;
  439.       a->lf = b;
  440.       b = a;
  441.       a = de;
  442.       de = deepcopy(b->rt);
  443.     }
  444.     else {                        /* a + b/c = (ac + b)/c */
  445.       b->kind = MUL;
  446.       de = b->lf;
  447.       b->lf = a;
  448.       a = b;
  449.       b = de;
  450.       de = deepcopy(a->rt);
  451.     }
  452.     nu = newoper(p->kind);
  453.     nu->lf = a;
  454.     nu->rt = b;
  455.     p->kind = DIV;
  456.     p->lf = nu;
  457.     p->rt = de;
  458.     ++r;
  459.   }
  460.   return r;
  461. }
  462.  
  463. /*-----------------------------------------------------------------
  464.    distribute2
  465.  
  466.    This applies the distributive laws
  467.    1. division over addition.
  468. */
  469. int distribute2(node *p) {
  470.   node *a,*b,*c;
  471.   int i,r=0;
  472.  
  473.   for (i=0; i<p->nump; ++i)
  474.     r+=distribute2(p->parm[i]);
  475.  
  476.   i=0;
  477.   c = p->parm[i];
  478.   if (p->kind==DIV && aop(c->kind)) {
  479.     b = deepcopy(p->parm[1-i]);
  480.     a = newoper(p->kind);
  481.     a->parm[i] = c->parm[1-i];
  482.     a->parm[1-i] = b;
  483.     c->parm[1-i] = p->parm[1-i];
  484.     p->parm[1-i] = a;
  485.     p->kind = c->kind;
  486.     c->kind = a->kind;
  487.     ++r;
  488.   }
  489.   return r;
  490. }
  491.  
  492. /*    This is used to control the direction of distribute */
  493. static int topdown = 0;
  494. /*-----------------------------------------------------------------
  495.    distribute
  496.  
  497.    This applies the distributive laws
  498.    1. multiplication over addition.
  499.    2. (ab)^2 = a^2*b^2
  500.    3. x^(a+b) = x^a*x^b
  501. */
  502. int distribute(node *p) {
  503.   node *a,*b,*c;
  504.   int i,r=0;
  505.  
  506.   if (!topdown)
  507.     for (i=0; i<p->nump; ++i)
  508.       r+=distribute(p->parm[i]);
  509.  
  510.   for (i=0; i<2; ++i) {
  511.     c = p->parm[i];
  512.     if (p->kind==MUL && aop(c->kind) ||
  513.         i && p->kind==EXP && aop(c->kind) ||
  514.         !i && p->kind==EXP && (c->kind==MUL || c->kind==DIV)) {
  515.       b = deepcopy(p->parm[1-i]);
  516.       a = newoper(p->kind);
  517.       a->parm[i] = c->parm[1-i];
  518.       a->parm[1-i] = b;
  519.       c->parm[1-i] = p->parm[1-i];
  520.       p->parm[1-i] = a;
  521.       p->kind = c->kind;
  522.       c->kind = a->kind;
  523.       if (i && a->kind==EXP) p->kind+=2;    /* change ADD to MUL */
  524.       ++r;
  525.       break;   /* just to be safe */
  526.     }
  527.   }
  528.  
  529.   if (topdown && !r)
  530.     for (i=0; i<p->nump; ++i)
  531.       r+=distribute(p->parm[i]);
  532.  
  533.   return r;
  534. }
  535.  
  536. /*-----------------------------------------------------------------
  537.    distribute_c
  538.    this applies the distributive law over multiplication which is
  539.    not governed by another multiplication or equality.
  540. */
  541. int distribute_c(node *p) {
  542.   int r=0;
  543.  
  544.   topdown = 1;
  545.   if (p->kind==EQU) {
  546.     r+=distribute_c(p->lf);
  547.     r+=distribute_c(p->rt);
  548.   }
  549.   else if (p->kind==MUL) {
  550.     r+=distribute(p->rt);
  551.     r+=distribute_c(p->lf);
  552.   }
  553.   else r+=distribute(p);
  554.   topdown = 0;
  555.  
  556.   return r;
  557. }
  558.  
  559. /*--------------------------------------------------------------------
  560.    fixassoc
  561.      + * are left assoc, ^ is right assoc.
  562. */
  563. int fixassoc(node *p) {
  564.   int i,r=0;
  565.   node *a,*b;
  566.  
  567.   for (i=0; i<p->nump; ++i)
  568.     r+=fixassoc(p->parm[i]);
  569.  
  570.   a = p->lf;
  571.   b = p->rt;
  572.   switch (p->kind) {
  573.   case ADD:
  574.   case MUL:
  575.     if (b->kind==p->kind ||              /* a+(b+c) = a+b+c */
  576.         b->kind==p->kind+1) {            /* a+(b-c) = a+b-c */
  577.       swingb;
  578.       i=b->kind; b->kind=p->kind; p->kind=i;
  579.     }
  580.     else break;
  581.     ++r; break;
  582.   case SUB:
  583.   case DIV:
  584.     if (b->kind==p->kind) {             /* a-(b-c) = a-b+c */
  585.       swingb;
  586.       --p->kind;
  587.     }
  588.     else if (b->kind==p->kind-1) {      /* a-(b+c) = a-b-c */
  589.       swingb;
  590.       ++b->kind;
  591.     }
  592.     else break;
  593.     ++r; break;
  594.   case EXP:
  595.     if (a->kind==EXP) {                  /* (x^y)^z = x^(y*z) */
  596.       swinga;
  597.       a->kind = MUL;
  598.     }
  599.     else break;
  600.     ++r; break;
  601.   }
  602.   return r;
  603. }
  604.  
  605. /*-----------------------------------------------------------------
  606.    get term, return p without coefficient
  607. */
  608. node *get_term(node *p, int oper, double *r) {
  609.   *r = 1.0;              /* default coeff */
  610.   if (p->kind==oper && p->rt->kind==NUM) {
  611.     *r = p->rt->value;
  612.     return p->lf;
  613.   }
  614.   return p;
  615. }
  616. /*--------------------------------------------------------------------
  617.    Combine all terms with common base within a clag.
  618. */
  619. int combine(node *p) {
  620.   int i,oper,r=0;
  621.   node *a,*b;
  622.   node *t1,*t2;
  623.   double c1,c2;
  624.  
  625.   for (i=0; i<p->nump; ++i)
  626.     r+=combine(p->parm[i]);
  627.  
  628.   a = p->lf;
  629.   b = p->rt;
  630.   oper = p->kind;
  631.   if (oper!=ADD && oper!=MUL)          /* not a clag */
  632.     return r;
  633.   i = 1;
  634.   if (a->kind!=oper) {
  635.     a=p; i=0;
  636.   }
  637.  
  638.   t1 = get_term(a->parm[i],oper+2,&c1);
  639.   t2 = get_term(b,oper+2,&c2);
  640.   if (equal(t1,t2)) {            /* same base, combine terms */
  641.     if (a->parm[i]==t1)
  642.       freetree(t1);
  643.     else {                    /* free one of the expressions */
  644.       freetree(b);
  645.       b = a->parm[i];
  646.       t2 = t1;
  647.     }
  648.     c2 += c1;
  649.     if (b==t2) {             /* make room for coeff */
  650.       b = newoper(oper+2);
  651.       b->lf = t2;
  652.       b->rt = newnum(0);
  653.     }
  654.     b->rt->value = c2;
  655.     if (c2==1.0) {            /* remove unit coeff */
  656.       freenode(b->rt);
  657.       freenode(b);
  658.       b = t2;
  659.     }
  660.     if (i) {                  /* cleanup... 2 stage */
  661.       nodecpy(p,a);
  662.       freenode(a);
  663.       p->rt = b;
  664.     }
  665.     else {                    /* 1 stage */
  666.       nodecpy(p,b);
  667.       freenode(b);
  668.     }
  669.     ++r;
  670.   }
  671.   return r;
  672. }
  673.  
  674.  
  675.  
  676. /*-----------------------------------------------------------------
  677.    Simplify
  678.    This function converts a rational expression to normal form.
  679.    Some of the attributes of normal form: no negative exponents,
  680.    terms are sorted over mul and add.
  681.    The first movenums pushes all numbers to the left so that
  682.    when the sortnode runs it pushes them to the right and combines
  683.    them.
  684. */
  685. void simplify(node *p) { simplify2(p,0); }
  686.  
  687. void simplify2(node *p, int slow)
  688. {
  689.   while (calcnode(p,1))     if (slow) slow=debug(p);
  690.   while (fixassoc(p))       if (slow) slow=debug(p);
  691.   while (nosubdiv(p))       if (slow) slow=debug(p);
  692.   while (fixassoc(p))       if (slow) slow=debug(p);
  693.   while (movenums(p,0,MUL)) if (slow) slow=debug(p);
  694.   if (sortnode(p,MUL) && slow) slow=debug(p);
  695.   while (movenums(p,0,ADD)) if (slow) slow=debug(p);
  696.   if (sortnode(p,ADD) && slow) slow=debug(p);
  697.   while (combine(p))        if (slow) slow=debug(p);
  698.   while (calcnode(p,1))     if (slow) slow=debug(p);
  699.   if (sortnode(p,MUL) && slow) slow=debug(p);
  700.   if (sortnode(p,ADD) && slow) slow=debug(p);
  701.   while (bisect(p))         if (slow) slow=debug(p);
  702.   while (calcnode(p,1))     if (slow) slow=debug(p);
  703.   while (movenums(p,0,MUL)) if (slow) slow=debug(p);
  704. }
  705.  
  706. /*-----------------------------------------------------------------
  707.    substitute -- prereq tgt MUST be an EQU.
  708. */
  709. void substitution(node *p) {
  710.   int i;
  711.   if (equal(p,tgt->lf))
  712.     movenode(p,deepcopy(tgt->rt));
  713.   else for (i=0; i<p->nump; ++i)
  714.     substitution(p->parm[i]);
  715. }
  716.  
  717.  
  718. /*-----------------------------------------------------------------
  719.    insert key
  720.    has special handling for equations...
  721. */
  722. void insertkey(int opr) {
  723.   node *tmp,*t1,*t2;
  724.  
  725.   if (!src || !tgt) return;
  726.   if (tgt->kind==EQU) {
  727.     t1 = tgt->lf; t2 = tgt->rt;
  728.   } else t1 = t2 = tgt;
  729.  
  730.   if (src->kind==EQU) {
  731.     tmp = src->lf;
  732.     src->lf = newoper(opr);
  733.     src->lf->lf = tmp;
  734.     src->lf->rt = deepcopy(t1);
  735.     tmp = src->rt;
  736.     src->rt = newoper(opr);
  737.     src->rt->lf = tmp;
  738.     src->rt->rt = deepcopy(t2);
  739.     src->rt->lf = tmp;
  740.   }
  741.   else if (opr==EXP) {
  742.     tmp = newnode();
  743.     nodecpy(tmp,src);
  744.     src->kind = EXP;
  745.     src->nump = 2;
  746.     src->lf = cons(tmp,EXP,deepcopy(t1));
  747.     src->rt = cons(newnum(1),DIV,deepcopy(t2));
  748.     src = src->lf;
  749.   }
  750.   else {
  751.     tmp = newnode();
  752.     nodecpy(tmp,src);
  753.     src->kind = opr ^ 1;
  754.     src->nump = 2;
  755.     src->lf = newoper(opr);
  756.     src->rt = deepcopy(t2);
  757.     src = src->lf;
  758.     src->lf = tmp;
  759.     src->rt = deepcopy(t1);
  760.   }
  761. }
  762.  
  763. /*-----------------------------------------------------------------
  764.    invert left side of equality, i=0 or 1
  765. */
  766. void cross_eq(node *p, int i) {
  767.   node *t;
  768.   if (p->kind != EQU) return;
  769.   t = p->lf;
  770.   if (t->kind==ADD ||
  771.       t->kind==SUB ||
  772.       t->kind==MUL ||
  773.       t->kind==DIV) {
  774.     p->lf = t->parm[i];
  775.     t->parm[i] = p->rt;
  776.     p->rt = t;
  777.     if (i && !(t->kind&1)) {   /* swap lf and rt */
  778.       t = t->rt;
  779.       p->rt->rt = p->rt->lf;
  780.       p->rt->lf = t;
  781.       t = p->rt;
  782.     }
  783.     if (!i || !(t->kind&1)) {   /* toggle operator */
  784.       t->kind ^= 1;
  785.     }
  786.   }
  787. }